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A statistical-mechanical investigation is performed on Rayleigh-Benard convection of a dilute 
classical gas starting from the Boltzmann equation. We first present a microscopic derivation of basic 
hydrodynamic equations and an expression of entropy appropriate for the convection. This includes 
an alternative justification for the Oberbeck-Boussinesq approximation. We then calculate entropy 
change through the convective transition choosing mechanical quantities as independent variables. 
Above the critical Rayleigh number, the system is found to evolve from the heat-conducting uniform 
state towards the convective roll state with monotonic increase of entropy on the average. Thus, the 
principle of maximum entropy proposed for nonequilibrium steady states in a preceding paper [T. 
Kita: J. Phys. Soc. Jpn. 75 (2006) 114005] is indeed obeyed in this prototype example. The principle 
also provides a natural explanation for the enhancement of the Nusselt number in convection. 



I. INTRODUCTION 

In a preceding paper^ we have proposed a principle 
of maximum entropy for nonequilibrium steady states: 
The state which is realized most probably among possible 
steady states without time evolution is the one that makes 
entropy maximum as a function of mechanical variables. 
We here apply it to Rayleigh-Benard convection of a di- 
lute classical gas to confirm its validity. 

Rayleigh-Benard convection is a prototype of nonequi- 
librium steady states with pattern formation, and exten- 
sive studies have been carried out to clarify it 
However, no calculation seems to have been performed 
on entropy change through the nonequilibrium "phase 
transition," despite the fact that entropy is the key 
concept in equilibrium thermodynamics and statistical 
mechanics. There may be at least four reasons for 
it. First, there seems to have been no established ex- 
pression of nonequilibrium entropy. Second, the stan- 
dard starting point to describe Rayleigh-Benard convec- 
tion is a set of deterministic equations for the particle, 
momentum and energy flows, with all the thermody- 
namic effects pushed into phenomenological parameters 
of the equationsi 2 -^ Third, one additionally adopts the 
Oberbeck-Boussinesq approximation to the equations in 
the conventional treatment Despite many theoreti- 
cal efforts over a long period ; 11 ' 12 ! 13 ] 14 ? 15 a well-accepted 
systematic justification for it seems still absent, thereby 
preventing a quantitative estimation of entropy change. 
Fourth, there is ambiguity on what to choose as indepen- 
dent variables of entropy for open systems. 

In a preceding paper jl we have derived an expression 
of nonequilibrium entropy together with the evolution 
equations for interacting bosons/fermions. We here ap- 
ply them to a classical gas of the dilute high-temperature 
limit where the evolution equations reduce to the Boltz- 
mann equation. We carry out a microscopic derivation of 
the hydrodynamic equations for the particle, momentum 
and energy densities (i.e., the basic conservation laws) 
from the Boltzmann equation. We then provide a system- 
atic justification for the Oberbeck-Boussinesq approxi- 
mation to describe the convection. With these prelim- 



inaries, we perform a statistical-mechanical calculation 
of entropy change through the convective transition by 
choosing mechanical quantities as independent variables. 
It is worth pointing out that classical gases have been 
used extensively for detailed experiments on Rayleigh- 
Benard convection over the last two decades£&2 Thus, 
quantitative comparisons between theory and experiment 
are possible here. 

This paper is organized as follows. Section HTl derives 
(i) equations of motion for the particle, momentum and 
energy densities and (ii) an explicit expression for the dis- 
tribution function /, both starting from the Boltzmann 
equation. Section|TTT]reduces the equations of fJTT]in a way 
appropriate to treat Rayleigh-Benard convection of a di- 
lute classical gas. This includes a systematic justification 
of the Oberbeck-Boussinesq approximation and a deriva- 
tion of the expression of entropy for convection. Section 
II VI transforms the equations of ^Illl into those suitable for 
periodic structures with the stress-free boundaries. Sec- 
tion [V] presents numerical results obtained by solving the 
equations of ^IVI It is shown explicitly that the principle 
of maximum entropy is indeed obeyed by the convection. 
Concluding remarks are given in WII 



II. DISTRIBUTION FUNCTION, 
CONSERVATION LAWS AND ENTROPY 

A. The Boltzmann equation and entropy 

We shall consider a monatomic dilute classical gas un- 
der gravity. This system may be described by the Boltz- 



mann equation 
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Here f = f(p, r, t) is the distribution function, t the time, 
p the momentum, m the mass, r the space coordinate, 
and g the acceleration of gravity. With a unified de- 
scription of classical and quantum statistical mechanics 
in mind, we choose the normalization of / such that it is 
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dimensionless and varies between and 1 for ferrnions. 
The collision integral C is given explicitly by 



C(p,r,t) 
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where V q is Fourier transform of the interaction potential 
and E p =p 2 /2m. 

Equation (TT]) also results from eq. (63) of ref. [l] as fol- 
lows: (i) approximate the spectral function as A(pe, rt) = 
2it8{e — Ep — mgz); (ii) substitute the second-order self- 
energy of eq. (66) into the collision integral of eq. (64); 
(iii) take the high-temperature limit; and (iv) integrate 
eq. (63) over e to obtain an equation for 



f(p,r,t) 



de 

—A(pe,rt)(j){pe,rt). (3) 



This whole procedure amounts to treating the interaction 
potential only as the source of dissipation in the dilute 
high-temperature limit, thus neglecting completely its in- 
fluence on the density of states, i.e., on the real part of 
the self-energy. 

With a change of variables p\ — p + q, p' = P — q'/2 
and p[ = P + q'/2 in eq. ([5]), the collision integral is 
transformed into 
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j [f p +(q-q')/2fp+(q+q')/2~ fpfp+q] , 
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where da — dfl q > § dq 1 8{q' —q) \m\V( q > _ q y 2 \l^h 21 \ is the 
differential cross section in the center-of-mass coordinate 
of the scattering with dfl q i denoting the infinitesimal 
solid angle. We shall use the contact interaction with 
no q dependence in V q , where da reduces to 



da = 



2 dn q , Jdq'S(q'-q), 



(5) 



with a = m\V\/4irh 2 denoting the scattering length. Now, 
the differential cross section acquires the form of the two- 
particle collision between hard-sphere particles with ra- 
dius a~2- 

Entropy per unit volume is given in terms of / by 
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with V denoting the volume. This expression also results 
from eq. (69) of ref. [l] for entropy density by: (i) adopting 
the quasiparticle approximation A(pe, rt) = 2it6(e— E p — 
mgz); (ii) defining / by eq. ([3]) above; (iii) taking the 
high-temperature limit; and (iv) performing integration 
over r. The term —1 in the integrand of eq. ([6]) is absent 
in the Boltzmann ii-functior>i£ but naturally results from 
the procedure (iii) above. Indeed, eq. ^ reproduces the 
correct expression of entropy in equilibrium. 19 



B. Conservation laws 

We next consider conservation laws which originate 
from the Boltzmann equation. The basic physical quan- 
tities relevant to them are the density n(r,t), the velocity 
v(r, i), the temperature T(r, t), the momentum flux den- 
sity tensor n(r, t) in the reference frame moving with the 
local velocity v, and the heat flux density Jq(t, t). They 
are defined by 
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with p = p — mv. The expressions (I7al) - (|7e)) also result 
from eqs. (C-la), (C-lb), (C-20a), (C-14) and (C-20b) of 
ref. Ill, respectively, by noting eq. ((3]), neglecting the in- 
teraction terms, and identifying the local internal-energy 
density £ as £ — \nk B T. 

To obtain the number, momentum and energy conser- 
vation laws, let us multiply eq. (fTJ) by 1, p and p 2 /2m, 
respectively, and perform integration over p. The con- 
tribution from the collision integral ^ vanishes in all 
the three cases due to the particle, momentum and en- 
ergy conservations through the collision] 16 i 17 The result- 
ing hydrodynamic equations can be written in terms of 
the quantities of eq. ([7]) as 
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where e z is the unit vector along the z axis and A : B_ 
denotes the tensor product: A:B_ = J^ij AijBji- These 
equations are identical in form with eqs. (C-2), (C-9) and 
(C-21) of ref. [l], respectively, with U = mgz and £ = 
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C. The Enskog series 

We now reduce the whole procedures of solving eq. ([1]) 
for (p, r,t) to those of solving eq. © for (r,t). We adopt 
the well-known Enskog metho d 16 ' 17 for this purpose, i.e., 
the expansion from the local equilibrium. We here de- 
scribe the transformation to the extent necessary for a 
later application to Rayleigh-Benard convection. 

Let us expand the distribution function formally as 



with P = nk B T the pressure and 1 the unit tensor. The 
left-hand side of eq. (fl2|) is thereby transformed into 
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with fc a dimensionless quantity defined by 



f(p,r,t) = f^(p,r,t) l + ^W(p,r,i) 



(9) 



where /( cq ) is the local-equilibrium distribution given ex- 
plicitly by 
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with p = p — mv. This /( oq ) has been chosen so as 
to satisfy the two conditions! 16 ' 17 (i) the local equilib- 
rium condition that the collision integral vanishes; (ii) 
eqs. (|7ap - ([7cl) by itself. It hence follows that the higher- 
order corrections < / 9 <J - ) {j = 1, 2, • • • ) in eq. ([9]) should obey 
the constraints: 



k = p/y / 2mk B T. 
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It is convenient to introduce the additional dimensionless 
quantities: 



2k B T 

and the mean-free path: 
/ 
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Using eqs. {TS} and (T7]) and noting eqs. ©, flD]). (H 
and (|15J) . we can transform eq. (|12p into the dimension- 
less form: 



(27Tft) 3 ^ J 
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Note that we have incorporated five space-time depen- 
dent parameters in f^ cq \ i.e., n, v and T, which can 
be determined completely with eqs. (|5a|) - ([5cl) of conser- 
vation laws. The remaining task here is to express the 
extra quantities II and jq in eq. ([8]) as functionals of n, 
v and T. 

Let us substitute eq. j9|) into eq. (P) , regard the space- 
time differential operators on the left-hand side as first- 
order quantities, and make use of the fact that the colli- 
sion integral ([2]) vanishes for /( cq ) . We thereby arrive at 
the first-order equation: 



9/(eq) p 3j(eq) 



at 



+ P*L f («Q =C (l) (12) 

or k B l 



2^ f&q fdV 8(q'-q) c _y_ Ch+Q? 
(i) , (1) _ (i) 

^fc+(q-q')/2 + ^fe+(q+q')/2 ^ ^+9 

2( fcfc-yl) :Vt> + fp-^fc-VlnT 



(18) 



where V = <9/<9r, and we have redefined tp^ as a function 
of k = p j \J2mk B T . Similarly, eq. ([TT]) now reads 



d^e-^fc™^^ = (n = 0,l,2). (19) 



Equation (fl8|) with subsidiary condition (I19|) forms a lin- 
ear integral equation for ip. . 

The right-hand side of eq. (fT8j) suggests that we may 
seek the solution in the form : 16 ' 17 



where is obtained from eq. ((H) as 
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With eq. (UHl) , the derivatives of / (oq) in eq. (H]) are 
transformed into those of n, v and T. We then remove the 
time derivatives by using eq. ([8]) in the local-equilibrium 
approximation, where 



n (oq) =pi, 



jj q) =°- 
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where ^4 5 / 2 and A 3 / 2 are two unknown functions; the use 
of fractions a = 5/2 and 3/2 to distinguish them will be 
rationalized shortly. Substituting eq. (|20[) into it, we can 
transform eq. (fl~8|) into separate equations for A a as 



2V2 f^q /dV %^9) e -P-(fc+q) 
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where tensor TZ% and Tu are denned in Table [T] together 
with another tensor . Since the factor e is present 
on the right-hand side of eq. (|2"Tj) , we expand A a (e) fur- 
ther in the Sonine polynomials Sf(e) a o 16 i 17 



A"(e) 



1=0 



c2Sf{e). 



(22) 



Use of two different complete sets {Sf}t (a = 5/2, 3/2) is 
only for convenience to transform the right-hand side of 
cq. (|2ip into a vector with a single nonzero element. With 
eq. (|20|) and (j22|) and the orthogonality of S 1 " (e) , we find 
that the constraint (I19|) reduces to the single condition: 



„ 3 / 2 



= 0. 



(23) 



We hence remove the 1 = term of a = 3/2 from eq. (|22|) 
in the subsequent discussion. We now take the tensor 
(a = 5/2) or the vector (a = 3/2) product of eq. (|21~j) with 
5"(fc 2 )W£ : /47r and perform integration over k. Equation 
(I2T1) is thereby transformed into an algebraic equation for 
the expansion coefficients {c"}^ as 



E 



T 



(24) 



The quantities I" e ,(k, q, q) and I" {l (k, q, —q) are obtained 
from eq. ()27|1 by removing the integral over df2 q < /47r and 
setting q' — > q and q' — > — q in the integrand, respectively. 
The first few series of eq. (|2l)l) are easily calculated 

-1/4, 



-5/2 _ ,7-3/2 _ 



J ll — L ) 2 01 

3/2 



5/2 _ ^-3/2 



analytically as 7^ 

Tj 5 / 2 = 205/48 and r 2 J 2 /z = 45/16, in agreement with 
the values given below eq. (10.21,3) of Chapman and 
Cowling*^ The matrix element for a general ££' can eval- 
uated numerically. With T^, and TZ" thus obtained, eq. 
([24| is solved by cutting the infinite series at a finite value 
£ c , and £ c is increased subsequently to check the conver- 
gence. Table [TT] lists the values of c" thereby obtained. 



Those of Cq 2 and c] 1 z are about 2% larger in magnitude 
than the analytic ones 5y/ir/A and — Iby/n/W with £ c = 
and 1 for a = 5/2 and 3/2, respectively. This rapid con- 
vergence as a function of £ c was already pointed out by 
Chapman and CowlingJ^ 

Substituting eqs. ©, (TO]) and JTOD into eqs. (7d| and 
(Tfe)) . we arrive at the first-order contributions to the mo- 
mentum flux density tensor and the thermal flux density 
as 
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with (W,TZ) = W : K and W-TZ for a = 5/2 and 3/2, 
respectively. Also, T e ", is obtained with a change of vari- 
able k^k — q/2 in the k integral as 



dfce~ 2fe2 fc 2 / dqe~1 2 / 2 q 3 / dq'S(q'-q) 



rra 
1 W 
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x(W£_ q/2) W^ q , /2 ). (27) 
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respectively, where v and re are the kinematic viscosity 
and the thermal diffusivity; 2 ^ respectively, defined by 



1, 2k B T 5/2 
4 V th 



-\J==cT 

6 V to 



(29a) 
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These quantities clearly have the same dimension. Using 
them as well as the specific heat at constant pressure Cp 
and constant volume Cy, we can introduce an important 
dimensionless quantity Pr = (y j n)(Cp /Cy) called the 
Prandtl number. 16 Adopting Cp/Cv = 5/3 of the ideal 
monatomic gas, we find P?* = 0.66 from eq. (|29[) and Table 
im which is in excellent agreement with the value 0.67 
for Ar at T = 273K. 17 Table [Hi] lists values of relevant 
thermodynamic and transport coefficients around room 
temperature at 1 atm for Ne, Ar and air. 

Thus, we have successfully expressed II and jq in terms 
of n, v and T as eqs. (fllj). (|28]l and ((29j) within the first- 
order gradient expansion. Now, eq. ([8]) with eqs. (fl4|). 



TABLE I: Quantities appearing in eq. (J5TJ. Here Sg (e) = 1 
and S"(e) = l+a-e. 

a W£ 71% Tk 

5/2 fefe-(fc 2 /3)I 2S£(k 2 )W£ A a (k 2 )W£ 
3/2 k -ST(fc 2 )W£ A a (k 2 )W% 



TABLE II: Values of c a t obtained by solving eq. ([24 

_ ry. jy. „a ,,a 

a Cq C X C 2 C3 C 4 

5/2 2.2511 0.1390 0.0233 0.0058 0.0018 
3/2 -1.7036 -0.1626 -0.0371 -0.0117 
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(|28[) and (|29|) forms a closed set of equations for the five 
parameters n, v and T incorporated in f cq (p,r,t). Af- 
ter solving them, we can obtain the distribution function 
f(p,r,t) by eqs. ©, (JTJJJ, (J2DJ, © and Table M and 
subsequently calculate entropy by eq. ([6]). 



the entropy change through the convective transition on 
a firm ground. 



Introduction of dimensionless units 



III. APPLICATION TO RAYLEIGH-BENARD 
CONVECTION 

We now apply the equations of EjTT] for n, v and 
T to Rayleigh-Benard convection of a dilute classical 
monatomic gas confined in the region —d/2 < z < d/2 
and —L/2<x,y<L/2. The gas is heated from below so 
that 

T(x,y,z = ±d/2) =T t AT/2, AT > . (30) 

The thickness d and the lateral width L are chosen as 
I <C d <C L. It hence follows that (i) there are enough 
collisions along z and (ii) any effects from the side walls 
may be neglected. We eventually impose the periodic 
boundary condition in the xy plane. 

We study this system by fixing the total particle num- 
ber, total energy, and total heat flux through z= —d/2. 
This is equivalent to choosing the average particle den- 
sity n, the average energy density £ , and the average 
heat flux density j'q at z = —d/2 as independent vari- 
ables; hence To =To(n,£, j'q) and AT = AT(n, £ ,j<g) in 
eq. (|3"D|) . The latter two conditions also imply, due to 
the energy conservation law, that there is average energy 
flux density jq through any cross section perpendicular 
to z. The fact justifies our choice of jQ as an independent 
variable to specify the system. It should be noted that 
this energy flow in the container may be due partly to a 
macroscopic motion of the gas. 

The standard theoretical treatment of Rayleigh- 
Benard convection starts from introducing the Oberbeck- 
Boussinesq approximation to the equations for n, v and 
T ^ However, this approximation seems not to have been 
justified in a widely accepted way. We here develop a 
systematic approximation scheme for the equations in fJII] 
appropriate to treat Rayleigh-Benard convection, which 
will be shown to yield the equations with the Oberbeck- 
Boussinesq approximation as the lowest-order approxi- 
mation. This consideration also enables us to estimate 



TABLE III: Values of relevant thermodynamic and transport 
coefficients under atmospheric pressure given in CGS units 
\a = V- Y (dV/dT)\. 



Ne (0°C) 



mn 0.900x10 
q 3.66x10" 
v 33.0x10" 

K 



83.1x10" 



Ar (0°C) 

I. 78x10 
3.67x10 

II. 8x10 
29.2x10 



air (0°C) 
1.29xl0~ 3 
3.67xl0~ 3 
13.2 xl0~ 2 
25.9 xlO -2 



We first introduce a characteristic temperature T de- 
fined by 



T = 2£/Mk B 



(31) 



We then adopt the units where the length, velocity and 
energy are measured by d, y/ k B T/m and k B T, respec- 
tively. Accordingly, we carry out a change of variables as 



t = d 



k B T 



and 
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dr'. 
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Let us substitute eqs. HU), © and (JUJ) into eq. © and 
subsequently perform the above change of variables. We 
thereby obtain the dimensionless conservation laws: 



2£ + V'.(nV) = 0, 



(33a) 
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with P' = n'T' and 
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An important dimensionless quantity of the system is the 
Rayleigh number R defined by 



R = 



_ U' g AV _ gf^ATd 3 



VK 



(35) 



where T _1 appears as the thermal expansion coefficient 
a of the ideal gas. 

The above equations will be solved by fixing n, £ = 
inknT/2 and jq, as already mentioned. These condi- 
tions are expressed in the dimensionless form as 



± |nV)dV = n'. 
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where P' = n'T', and integrations extend over —L'/2 < 
x',y' < L'/2 and -1/2 < z' < 1/2. Equation ([36b"f has 
been obtained by integration of (p 2 /2m + mgz)f over r 
and p with eq. l[7]). whereas eq. (|36c|) originates from eq. 
(l28bl) . 

To make an order-of-magnitude estimate for the pa- 
rameters in eqs. ([34)) and (|35|) . consider Ar of 273K at 
1 atm confined in a horizontal space of d cm with the 
temperature difference ATK. Using Table HT1 we then 
obtain the numbers: 



v' = 4.95 x l(T 6 /(i, k 1 = 1.22 x lCT 5 /d, 
= 1.73 x lCT 6 d, AT' = 3.67 x 1(T 3 AT, 



and 



R = 1.04 x 10 2 d s AT. 



(37a) 



(37b) 



The critical Rayleigh number R c for the convective tran- 
sition is of the order 10 3 ; 2 - which is realized for 2cm 
and AT~ IK. We now observe that the dimensionless pa- 
rameters have the following orders of magnitude in terms 
of (5=1CT 3 : 



AT' 



R 



(38) 



Thus, Rayleigh-Benard convection is a phenomenon 
where two orders of magnitude (i.e., 6 and 5 2 ) are rele- 
vant. From now on we shall drop primes in every quantity 
of eqs. J33J and (55)) . 



B. Omission of the number conservation law 

Let us write eq. (|3"3"|) in terms of j = nv instead of v. 
It follows from the vector analysis that the vector field 
j can be written generally as j = Vf+Vxj4, where $ 
and A corresponds to the scalar and vector potentials of 
the electromagnetic fields, respectively. We then focus in 
the following only on those phenomena where the current 
density satisfies V<& = 0, i.e., 



Vj = 0. 

This implies that we may drop eq. Q33aft from eq. 
treat only eqs. (|33b[) and 



C. Expansion in 8 



(39) 



to 



Equation (|3"5|) suggests that we may solve eqs. (|33bp 
and (|33cp in powers of S. Noting eqs. (|5oa)> and ((36b)) . 



we first expand n and T as 



n = n 1 



E 



T = 1 



(40a) 



e=i 



With eqs. ([291) . (I34p and ([38)) . we next expand dimension- 
less parameters v = {l/4)V2Tc 5 /2 , k = -(5l/6)V2Tcl /2 



and U g as 



« = £«W C/ S = lf>, (40b) 



£=2 



£=2 



where z/< 2 ) = (//4) % /2 C g /2 and k^ 2 ) = -{U /Q)y/2c\ /2 are 
constants with l — l(n). It also follows from AT~<5 and 
eq. (|36cl) that 



-.(3) 



(40c) 



It remains to attach orders of magnitude to the differen- 
tial operators and j = nv. In this context, we notice that 
the Oberbeck-Boussinesq approximation yields a critical 
Rayleigh number R c which is in good quantitative agree- 
ment with experiment.— The fact tells us that the pro- 
cedure to attach the orders should be carried out so as 
to reproduce R c of the Oberbeck-Boussinesq approxima- 
tion. The requirement yields 



^ 3 ' at 



0{5 



1.5\ 



V = O(S~ - 2b ). 



(40d) 

See eq. ([55]) below and the subsequent comments for de- 
tails. The above power-counting scheme will be shown 
to provide not only a justification of the Oberbeck- 
Boussinesq approximation but also a systematic treat- 
ment to go beyond it. 

Let us substitute eq. (@DJ to eqs. ((33b)) and ([36b)) . The 
contributions of O(S) in these equations read Vf' 1 ' =0 
and / pWd 3 r = 0, respectively, with pW = n(n^+T^). 
We hence conclude PW=0, i.e.. 



>(i) 



-T (1) 



(41) 



It also follows from eqs. (|36a|) and (|36c|) with eqs. ([40)) 
and (JUJ) that TW should obey 



T (1) d 3 r = 0, 



1 rL/2 g T (l) 

L 2 J- L /2 dX J~L/2 dy 9z 



= -1/2 



(42a) 



(42b) 



Equation ([4"2"]) is still not sufficient to determine It 
turns out below that the required equation results from 
the 0(5 3 ) and 0(6 2 5 ) contributions of eqs. (|33b|) and 
([33c[) . respectively. 

Next, collecting terms of 0(S 2 ) in eqs. (|33b[) and (|36b|) 
yield VP^ - 



-nU^ 2) e z 



and J(|p( 2 )+nC/^ J z)d 3 r 
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0, respectively. Hence P^ — —nUgZ. Noting — 

(l) T (l) +T (2) ) and 

using eq. (|4"Tj) . we obtain 
= (T«) 2 -Uf>z. (43) 

It follows from eqs. (|36ap and (I36c|) with eqs. ([3D]) and 
(|43| that T( 2 ) should obey 

[( T a))2_ T {2)] d 3 r = 0) (44a) 

= 0. 

(44b) 

In deriving eq. (|44b[) . use has been made of {nn)^ = 
ntt^'T^ 1 ' /2 which results from koc ZT 1 / 2 and 2 cxn , see 
eqs. (JTrJ) and (j29b|) . Equation (j44|) forms constraints on 
the higher-order contribution T^ 2 \ which will be irrele- 
vant in the present study, however. 

Finally, we collect terms of 0(5 3 ) in eq. (|33b[) to obtain 



,L/2 




/ dx 




l-L/2 J 


-L/2 V 



aT (2) ji(l) 5T (1) 



C*2 



z = -l/2 



0* 



:-(i-5) 



n«[/( 2 > 



VP< 3 > 



?; 
0. 



(45) 



where we have used eq. (|39p . We further operate VxVx 
to the above equation and substitute eq. (l4Tj) . This yields 



d_ 

dt 



V 2~-(i.5) + v x V x (j^ 1 - 5 ) • Vj' (1 - 5) ) 



W 2) (V 2 ) 2 ^ 1 ' 5 ) + C/f ) (e, V 2 - ■ VV)T« = . (46a) 
On the other hand, terms of 0(<5 2 5 ) in eq. (I33c|l lead to 
9T« 



■ j'Ci-B) . VT« - K < 2 )V 2 TW =0. 



(46b) 



Equation (|46p forms a set of coupled differential equa- 
tions for and j^" 5 \ which should be solved with eq. 
(|42p . It is almost identical in form with that derived with 
the Oberbeck-Boussinesq approximation, predicting the 
same critical Rayleigh number R c as will be shown below. 
The whole considerations on Raylcigh-Benard convection 
presented in the following will be based on eq. (|4*6]> with 
eq. (Hg). 

Two comments are in order before closing the subsec- 
tion. First, if we apply the procedure of deriving eq. 
(1161 t0 the °( s4 ) and 0(S 3 5 ) contributions of eqs. (|33b)) 
and (|33c|) . respectively, we obtain coupled equations for 
the next-order quantities and J^ 2 ' 5 ', which should 
be solved with eq. (|44| . Thus, we can treat higher- 
order contributions systematically in the present expan- 
sion scheme. Second, cq. (l4"5j) may be regarded as the 
equation to determine P^ for given and jt 1 - 5 ). It 
yields a relation between and iv- 3 ', which in turn 
leads to the constraint for upon substitution into 
eqs. (|36a[) and (|36c|) . On the other hand, the equation 
for originates from the 0(S 5 ) and 0(5 4 5 ) contribu- 
tions of eqs. (|33b[) and (|33cp . respectively. Now, one may 
understand the hierarchy of the approximation clearly. 



D. Expression of entropy 

We now write down the expression of entropy in pow- 
ers of 5. Entropy of the system can be calculated by 
eq. ([5]), where the distribution function / is given by 
eq. © with eqs. (fTU|) and pp|) . Hereafter we shall drop 
the superscript in ip^\ which specifies the order in the 
gradient expansion, to remove possible confusion with 
the expansion of eq. (|4TJ)l . Thus, / is now expressed as 
f = f^(l+<p). 

We first focus on ip and write eq. P0|) in the present 
units with noting eq. (|16p . We then realize that (p is pro- 
portional to (VT or ZVu, which are quantities of (3(<5 3 ) 
and 0(S 35 ) in the expansion scheme of eq. (|4"0|) , respec- 
tively. It hence follows that there is no contribution of 
0(5 2 ) from tp. In contrast, /( eq ) yields terms of 0(5 2 ), 
as seen below. Thus, we only need to consider /( cq ). 

Let us write /( eq ) of eq. (fTU|) in the present units, sub- 
stitute eq. (I40p into it, and expand the resulting expres- 
sion in powers of S. We also drop terms connected with 
j (i.e., v) which have vanishing contribution to S within 
0(5 2 ) after the momentum integration in eq. ([6]). We 
thereby obtain the relevant expansion: 

/ M = f|f^-'(i+|> ( '>) 



where e=p 2 /2, and and are defined by 
3 



w (1) (e) = e 



(48) 



Let us substitute eq. (147)) into eq. (HJ) and carry out in- 
tegration over p. The contribution of 0(1) is easily ob- 
tained as (fee = 1) 



In 



(2?r 



,3/2 



(2nh) 3 n 



(49) 



which is just the equilibrium expression^ for density n 
and temperature T in the conventional units, as it should. 
Next, we find =0 due to eqs. (gl]) and (J42a) . Thus, 
the contribution characteristic of heat conduction starts 
from the second order. A straightforward calculation 
yields 



s m 



(50) 



Equation ([50|) is the basic starting point to calculate the 
entropy change through the convective transition. Note 
that we have fixed jq in the present consideration, i.e., 
the initial temperature slope as seen from eq. (|42b|) . 

At this state, it may be worthwhile to present a qualita- 
tive argument on entropy of Rayleigh-Benard convection. 
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With the initial temperature slope fixed as eq. (|42b|) . eq. 
([50)1 tells us that entropy will be larger as the temper- 
ature profile becomes more uniform between z — ±1/2. 
The conducting state with v = has the linear tempera- 
ture profile, as shown shortly below in eq. (|52p . Thus, any 
increase of entropy over this conducting state is brought 
about by reducing the temperature difference between 
2 = ±1/2. Such a state necessarily accompanies a tem- 
perature variation which is weaker around z = than 
near the boundaries 2 = ±1/2. This temperature profile 
is indeed an essential feature of Rayleigh-Benard convec- 
tion which shows up as an increase of the Nusselt number 
(i.e., the efficiency of heat transport) under fixed temper- 
ature difference. 2 - 7 Combining eqs. (|42bj) and (f50|) with 
the experimental observation on the Nusselt number, we 
thereby conclude without any detailed calculations that 
entropy of Rayleigh-Benard convection should be larger 
than entropy of the conducting state in the present condi- 
tions with j'q =const. Thus, Rayleigh-Benard convection 
is expected to satisfy the principle of maximum entropy 
given at the beginning of the paper. We shall confirm 
this fact below through detailed numerical studies. 

IV. PERIODIC SOLUTION WITH 
STRESS-FREE BOUNDARIES 

Equation (jlS]) with eq. (f4"2"| forms a set of simulta- 
neous equations for T^ 1 ' and j^' 5 \ which should be 
supplemented by the boundary condition on For 
simplicity, we here adopt the assumption of stress-free 
boundaries £ 

= |^ 5) = at z = ±\. (51) 

However, qualitative features of the convective solutions 
will be universal among the present and more realis- 
tic/complicated boundary conditions; see the argument 
at the end of the preceding section. 

We first discuss the heat-conducting solution of eq. 
and its instability towards convection. We then trans- 
form eq. with eqs. (f4"2")) and (fSTj) in a form suitable 
to obtain periodic convective structures. 

A. Conducting solution 

Let us consider the conducting solution of eq. (j46|) 
where jt 1 - 5 ) =0 with uniformity in the xy plane. Equa- 
tion (|46|) then reduces to d 2 T*- 1 Ydz 2 = 0, which is solved 
with eq. f42]) as 

T« = — AT hc z , AT hc ee ^g^y . (52) 

Substituting this expression into eq. (|50p . we obtain en- 
tropy of the conducting state measured from 

5 (o) 

as 

^-T^c) 2 . (53) 



B. Instability of the conducting state 

We next check stability of the conducting solution by 
adding a small perturbation given by 

j( x - 5 )(r,t) = e xt e ik ^ r (6f smk z ( + 6jlcosk z () , (54a) 

T« (r,t) = —AT hc 2 + STe xt e lk± ' r sin k z ( , (54b) 

where £ ee z + 1/2, k z = £ 3 ir (£ 3 — 1,2, •••) from eq. 
(f5"Tj) . and 5j± denotes a vector in the xy plane. Let 
us substitute eq. (f54|) into eq. (|46|) and linearize it with 
respect to the perturbation. This leads to 

{\+v^k 2 )k 2 5f - Uf ) k\8Te z = , (55a) 
{\ + v {2) k 2 )k 2 5jl - iU^ 2) k ± k z ST = , (55b) 

(A + K (2) fc 2 )(5T - AT hc 5f z = . (55c) 

The components Sj^ and Sj^ in the xy plane are ob- 
tained from eqs. (|55a|) and (|55bj) as 

« = "• « = (irw^ < 56 > 

In contrast, the z component of eq. (|55ap and eq. (|55c|) 
form linear homogeneous equations for 8j s z and ST. The 
requirement that they have a non-trivial solution yields 

A 2 +(z/( 2 )+ K ( 2 ))fc 2 A+^ 2 )K( 2 )(fc 4 -i?(- 1 )fc 2 _/fc 2 ) = 0, (57) 

with R}- -1 ' the Rayleigh number defined by eq. (|35|) with 

U' g -> U ( g 2 \ AT -> AT hc , v' -> z/ 2 ) and k' -> . The 
conducting solution becomes unstable when eq. (|57p has 
a positive solution, i.e., 

_ £/f AT hc (fcj + fc 2 ) 3 27tt 4 
R = K ( 2 M 2 ) ~ k\ (58) 

Thus, we have obtained the value R c = 277r 4 /4 for 
the critical Rayleigh number— which corresponds to 
(k ± ,k z ) = (n/V2,ir). 

Besides reproducing the established results^ the above 
consideration may also be important in the following re- 
spects. First, we require that: (i) terms of the z com- 
ponent of eq. (|55a[) all have the same order in S with 
ST = O(S); (ii) the same be true for terms of eq. (|55c|) . 
This leads to the attachment of the order-of-magnitudc: 
5j = 0(S 1 - 5 ), fc = 0(<5-°' 25 ) and A = 0(<5 15 ), thereby justi- 
fying eq. (|40d[) . The conclusion k = O(S~ 025 ) also results 
from k = ^3/2tt ~ 4 for the critical Rayleigh number. 
Second, eq. ([55]) removes the ambiguity in v, k and a to 
estimate the critical Rayleigh number R c . Specifically, 
we should use the mean values over — l/2<z<l/2 for 
a detailed comparison of R c between theory and experi- 
ment. 
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C. Convective solution 

We now focus on the convective solution of eq. ([46]) 
with periodic structures. Let us introduce basic vectors 
as 



1 — cos k z 



(64b) 



a% = (ai x ,ai v ,0), a 2 = (0,a 2) 0), a 3 



(59) 



We consider the region in the xy plane spanned by A/iai 
and Mian with Mj (j = 1,2) a huge integer, and im- 
pose the periodic boundary condition. The wave vector 
k is then defined in terms of the reciprocal lattice vectors 
bi = 27r(a 2 Xa3)/[(oiXa 2 )-03], b 2 = 2n(a 3 xa 1 )/[(a 1 x 
a 2 )-a3] and b 3 =7re z as 



where ATh c denotes temperature difference of the con- 
ducting state given explicitly in eq. (|52p . 

It follows from the energy conservation law that eq. 
(|42bj) should also hold at z — 1/2, which leads to an 
alternative expression for AT CV . Subtracting it from eq. 
(|64aj) yields the identity obeyed by Tfe x =o,fe^ 



^Tfe x =o,fe a fc z (l-cos k z ) = 0, 



(65) 



ft* 



k — ^ £jbj , 



(60) 



with £j denoting an integer. The analysis of ^IVBI sug- 
gests that the stable solution satisfies |6i| ~ \bi\^ix I \]~2. 

Equation (|5rJ|) tells us that Sj± and ST are out-of- 
phase. Also noting eqs. (|42[) and (f5T|) . we now write 
down the steady solution of eq. (|46[) in the form: 

J_L M= X! ^2j±ksin{k ± - r)cosk z ( , (61a) 
fcx(^o) fc* 

jf' 5) W= E E^ fcC0S ^' r ) sinfc ^' ( 61b ) 

fcj_(#0) fc, 



which is useful to check the accuracy of numerical calcu- 
lations. 

Let us substitute eq. (fBTj) into eq. (|4"6")l . A straightfor- 
ward calculation of using eq. ([63)) then leads to coupled 
algebraic equations for T k , j zk and j pk as 

K (2) fc 2 T fc - AT cv j zk 



k'k" 



— 5k'l+k' ± kx^K-K k " 



-5k>J_+k' ± k ± Sk<J+k' z k z + fl ^2^l f ^k'J_-k' ± kj_^k'J+k' z k z 



k± • j±k" k z j zk 
Sk ± o 



>k' +k' 



' k ± Sk 



f 1 ^^^^fe' x -fe^fex^i'-fci*:* + bk'l-k' ± k^k> z -k'< 



T {1) (r) =J2J2^ k cos ( k ±- ' r ) sink z ( - AT cv z - T x , ~^-k'lkJk'i+k' z k z 



= 0. 



(66a) 



k± fc. 



(61c) 

where £=2+1/2, and we have chosen jz" 5 ^ and as 
even functions in the xy plane without losing the gener- 
ality. The k summations in eq. (|6ip run over 



(62) 



-oo<£i<-l, l<l 2 <oo, 1<4<oo, 
0<h<oo, 0<£ 2 <oo, 1<£ 3 <oo, 



which covers all the independent basis functions. The 
condition (1391) is transformed into 



k-j±k + k z j zk = . (63a) 
Thus, j±k may be written generally as 

3±k = - p 3zk H j: jpk ■ (63b) 



+ 4 E 3 zk ' [ fc ' 3k"(Sk'l+k' ± k^k'J-kik z 
k'k" 

+5 k '^- k ' ±k± Sk'J+k' z k z + S k '^ +k ' ±k± Sk' z ' + k' z k z 
— ^k'j_-k' ± k±^k'J-k' z k z ~ ^k'j_-k'j_k±^k' z -k'Jk z 

+ (k± -j±k" ~ kg,jzi t ")(Sk' x -k'lkxfa'J-k'sk z 
+8 k >^- k±k± Sk> z -k'^kz + Sk^+k'ikx^K-k'jkz 

~^k' ± -k'lk±^k'J+k' z k z )] 



fc'fc" 



The constants AT CV and Ti in eq. (|61cP denote the tem- 
perature difference between z — ±1/2 and the average 
temperature shift, respectively. They are fixed so as to 
satisfy eq. (14"2")) as 



AT rv = ATI, 



E 



T k± =o,k z k z 



(64a) 



x (^k'J_+k' ± k ± ^k' z '-k' z k z ~ $k'l -k' ± k ± ^k'J+k' z k z ) 

+k ■ jk'k ■ jk"{Sk'i+k' ± k ± Sk'j+k' z k z 

— fik'l-k' ± k ± Sk'J-k z k z — §k'j_-k'lk±fik' z -k'jkj 

—(k± ■j±k> - k z j zk >)(k^ -j±k" - kzjzk") 

x (Sk' ± -k'^k ± S k 'J-k' z k z + &k'l-k' x k_i_5k z -k'£k z ) 

+k ■ j k '{k±_ ■ jj_ k » - k z j zk »)(S kl±+k ^ k± S k ' z -k' z 'k z 
-5k>, -k'lk±^k' z +k'jkj] = 0, (66b) 
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(2) .2~ , 1 ( e z x fej-) ' JXfe' r, ~ 

fe'fe" - 1 

+3k'l+k' ± k ± 3k' : ! +k'^k z — Sk'l-k' ± k ± fik' : ! ~k' z k z 
X OV^-fe'Xfex^fcJ.'-S^fc* + Sk'l-k' ± k ± fik<z-k'Jkz 

—fik'j_+k'lk ± $k' z -k%kz+fik' ± -k'lk x fiK'+K k ^] = ®> (66c) 



with ji^fc given by eq. (|63b|) . Finally, entropy in convec- 
tion is obtained from eqs. (|50"|) and (I61c| as 



9(2) _ 



+AT CV Tf. ±= o,k 



fe 

1+cos fc z 



(67) 



V. NUMERICAL RESULTS 

We now present numerical results on entropy of 
Rayleigh-Benard convection obtained by solving eqs. 
(|66a[) - (|66cl) . As a model system we consider Ar at 
T = 273K under atmospheric pressure and use the values 
((37)1 for the parameters in eqs. (|66a[) - (|66c|) . We also fix 
the heat flux density jq at z — —d/2 so that tempera- 
ture difference AX^c of the conducting state, expressed 
in terms of jq as eq. (f5"2j) . takes the value AT\ 1C = IK. 
The Rayleigh number is then controlled by chang- 

ing the thickness d. We here adopt the units described 
above eq. (pH?)) with fcs = 1- 

Recent experiments on Rayleigh-Benard convection 
have been performed most actively with compressed clas- 
sical gases^^ where controlled optical measurements of 
convective patterns are possible, ft has been pointed out 
that the basic n (oc Z _1 ) and T dependences of eq. (f29|) are 
well satisfied even for those gases under high pressures. 8 
Thus, the present consideration with eq. (|37p has direct 
relevance to those experiments. 



A. Numerical procedures 

To solve the coupled equations numerically, we first 
multiply eqs. p^ - plc]) by 10 8 , 10 10 and 10 10 , respec- 
tively, and rewrite them in terms of T' k = 10 3 Tfc and 
j' k = 10 5 j k to obtain equations of O(l). We next replace 
zeros of the right-hand sides by —dT k /dt', —dj' zk /dt' and 
—dj' k /dt', respectively. These time derivatives are what 
one would get by retaining time dependence in the expan- 
sion coefficients of eq. (foTj) . We then discretize the time 
derivatives as df k (t')/dt' « [f k (t' + At') - f k (t')]/At', 
for example. The summations over k are truncated by 
using a finite value £ c (> 5) in place of oo in eq. ([6^)1 . 



As for periodic structures, we investigate the three can- 
didates: the roll, the square lattice and the hexagonal 
lattice with \bi | = ~7r/\/2. With these preliminaries, 
we trace time evolution of the expansion coefficients un- 
til they all acquire constant values. Choosing At' < 0.005 
and £ c > 5 yields excellent convergence for the calcula- 
tions presented below. The initial state is chosen as the 
conducting state with small fluctuations T' k ~ 10~ 2 for 
the basic harmonics k. The constants AT CV and T\ have 
been updated at each time step by using eq. Also 
evaluated at each time step is entropy measured with 
respect to the heat-conducting state: 



(2) 

he 



(68) 



where S^J and Set' are given by eqs. (|53| and (|67f . re- 
spectively. We thereby trace time evolution of AS simul- 
taneously. The above procedure is carried out for each 
fixed periodic structure. 

We have studied the range: 1 < R^/ R c < 10. 
Although the region extends well beyond the Busse 
balloon 3,6 of stability for classical gases^* 8 -^ it will be 
worth clarifying the basic features of steady periodic so- 
lutions over a wide range of the Rayleigh number. 



? (2) 



B. Results 

Figure Q] plots AS as a function of t' for the Rayleigh 
number Ry^ 1 ' = 1.2R C which is slightly above the criti- 
cal value R c — 27tt 4 /A. The letters r, s and h denote (r) 
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FIG. 1: Time evolution of entropy measured with respect 
to the heat-conducting state for — 1.2R C . The letters 

r, s and h denote roll, square and hexagonal, respectively, 
distinguishing initial fluctuations around the heat-conducting 
solution; see text for details. The final state of t' > 80 is 
the roll convection, whereas the intermediate plateaus of s 
and h correspond to the square and hexagonal convections, 
respectively. 
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FIG. 2: Time evolution of entropy AS. The four curves cor- 
respond to the different Rayleigh numbers: i?' -1 - 1 = 1.2i? c , 
2.0.Rc, 5.0i? c and W.0R c . The initial state is the heat- 
conducting state with the fluctuation T'[1,0, 1] = 1.0 x 10~ 2 
and bi = 7r/v2, whereas all the final states are the roll con- 
vection. The broken line near the top indicates the upper 
bound of AS. 



roll, (s) square and (h) hexagonal, respectively, distin- 
guishing initial conditions. Writing T^=T'[£i, and 
introducing the angle 9 by 9 = cos -1 (£>i ^2), those initial 
conditions are given explicitly as follows: (r) T"[l, 0, 1] = 
1.00 x 10~ 2 with |&i|=tt/V2; (s) f '[1, 0, 1] = 1.01 x 10~ 2 
and T"[0, 1,1] = 0.99 x 10~ 2 with |6 X | = |6 2 | = ir/V2 
and 9 = tt/2; (h) T" [1,0,1] = f '[0,1,1] = 1.00 x 10~ 2 
and f '[1,1,1] = 1.01 x 10~ 2 with |bi| = |6 2 | = tt/\/2 
and 9 = 2ir/3. We observe clearly that entropy in- 
creases monotonically in all the three cases to reach a 
common final value, which is identified as the roll con- 
vection. Thus, the law of increase of entropy is well 
satisfied even in the open driven system of Rayleigh- 
Benard convection under the condition of fixed mechan- 
ical variables, i.e., particle number, volume, energy and 
energy flux. The intermediate plateaus seen in s and 
h correspond to the metastable square and hexagonal 
lattices, respectively, with (s) T"[l, 0, 1] ~ T'[Q, 1, 1] and 
(h) f'[l, 0, 1] ~ f'[0, 1, 1] ~ T"[l, 1, 1]. The escapes from 
those lattice structures are brought about by the tiny ini- 
tial asymmetry between (s) T'[l, 0, 1] and T'[0, 1, 1] and 
(h) f '[1,0,1] = T'[0,1,1] and f '[1,1,1], to form even- 
tually the roll convection of (s) T' [0,1,1] = and (h) 
f'[l,0, 1] = T' [0,1,1] = 0. Indeed, those escapes are 
absent if we choose (s) T"[1,0,1] = f '[0,1,1] and (h) 
f'[l,0,l]=f'[0, l,l]=f '[1,1,1] at t' = 0. We have also 
performed the same calculations for different initial val- 
ues of f'[l, 0, 1], f"[0, 1, 1], f'[l, 1, 1] = 10~ 6 ~10- 2 . The 
smaller initial fluctuations merely delay the first steep 
rise in AS with no observable quantitative changes there- 
after. 
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FIG. 3: Profile of the average temperature variation T(z) in 
the roll convection normalized by the temperature difference 
ATkc in the heat-conducting state. The Rayleigh numbers 
are = R c , 1.2R C , 2.0R C , 5.0R C and 10.0i? c from top to 

bottom on the left part. 



Thus, we have confirmed that the roll convection is sta- 
ble for the infinite horizontal area, which was predicted 
originally by Schluter, Lortz and Bussed for R ~ R c 
based on the linear stability analysis; see also refs.0 and 
7. It should be noted at the same time that the en- 
tropy differences among different structures are rather 
small. We hence expect that: (i) the order of the 
stability may easily be changed by finite-size effects, 
boundary conditions, etc.; (ii) initial conditions, fluc- 
tuations and defects play important roles in Rayleigh- 
Benard convection. These are indeed the features ob- 
served experimentally— 

Figure [5] plots time evolution of AS for four different 
Rayleigh numbers, all developing from the initial fluctu- 
ation T'[l, 0,1] = 1.00 x 10~ 2 with |6i|=7r/V2. Each fi- 
nal state is the roll convection, which is stabilized faster 
as becomes larger. The increase in AS is seen 

quite steep for R(~^ = 10R C followed by a small oscilla- 
tion. This oscillation in AS may be due either to (i) the 
fluctuations in particle number, momentum, energy and 
energy flux inherent to open systems or (ii) the initial 
correlations which causes the anti-kinetic evolution i 21 i 22 
Such fluctuations are also observed in a numerical study 
by Orban and Belleman o 21 i 22 for an isolated system and 
not in contradiction with the law of increase of entropy. 
We observe clearly that the principle of maximum en- 
tropy proposed at the beginning of which is relevant 
to the final steady state without time evolution, is in- 
deed satisfied here. The dotted line in Fig. [2]is the upper 
bound of entropy in convection, as may be realized from 
eq. (|50j) . As the Rayleigh number is increased further, 
entropy differences between different structures become 
smaller so that the system will eventually fall into the 
region where fluctuations dominate with no stable struc- 
ture. The turbulence observed in this region2i&£ may be 
connected with this instability. 
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FIG. 4: The length of the stable roll convection as a 
function of the normalized Rayleigh number ii' -1 ' /R c . Here 
|6i| = tt/V§ at i? ( - 1) /i?c = l 



example, the stable state for w- ^ / R c = 10 corresponds 
to |6i| =0.829, but entropy is increased by only 2x 10~ 9 
from the value AS" = 1.296 x 10~ 6 at \b 1 \=Tc/V2. Hence 
it is expected that (i) initial conditions, fluctuations and 
boundary conditions play important roles and (ii) it may 
take a long time, or even be impossible in certain cases, 
to reach the stable state. These are indeed in agreement 
with experimental observations^^ 

The whole our results have been obtained for the 
stress- free boundary condition (|51[) . However, the quali- 
tative results will be valid also for more realistic bound- 
ary conditions. Indeed, experiments on Rayleigh-Benard 
convection all exhibit the temperature profile with a steep 
change near the boundaries followed by moderate varia- 
tion in the bulk region as reflected in the enhancement 
of the Nusselt number^i And it is this feature in the 
present study which has caused increase of entropy in 
convection over that in the conducting state. 



Figure [3] shows profile of the average temperature vari- 
ation T(z) along z in the roll convection for five different 
Rayleigh numbers. Temperature has less variation as the 
Rayleigh number becomes larger, thereby increasing en- 
tropy of the system. Thus, we may attribute the forma- 
tion of convection to its efficiency for increasing entropy 
under fixed inflow of heat, i.e., the initial slope of temper- 
ature. Experiments on Rayleigh-Benard convection have 
naturally been carried out by fixing the temperature dif- 
ference rather than the inflow of heat, and formation of 
the convection has been discussed in terms of the change 
in the Nusselt number— (i.e., the efficiency of the heat 
transport) due to the increase of the initial temperature 
slope through the convective transition. Here we have 
seen that the same phenomenon can be explained with 
respect to the basic thermodynamic quantity of entropy. 
Thus, the principle of maximum entropy partly justifies 
the maximum heat transfer hypothesis by Malkus and 
Veronis^ 

All the above calculations have been carried out by 
fixing the periodic structure. They clearly show that the 
principle of maximum entropy is indeed obeyed. We now 
take the principle as granted and use it to determine the 
stable lattice structure. We carry this out within the roll 
convection by changing the value |&i|. 

Figure 2] plots |&i| as a function of normalized Rayleigh 
number i?( _1 )/i? c . As seen clearly, |&i| increases gradu- 
ally from the value n/y/2 at R^~^ /R c = 1. This tendency 
is in qualitative agreement with the experiment by Hu, 
Ecke and Ahlers^ using a circular cell of L/d^> 1; see 
also ref. @. It should be noted at the same time that the 
entropy difference around |&i| =7t/a/2 is quite small. For 



VI. CONCLUDING REMARKS 

The present study shows unambiguously that the prin- 
ciple of maximum entropy given at the beginning of Sjl]is 
indeed satisfied through the Rayleigh-Benard convective 
transition of a dilute classical gas. The result is encour- 
aging for the principle as a general rule to determine the 
stability of nonequilibrium steady states. We need to in- 
vestigate other open systems, as well as Rayleigh-Benard 
convection without fixing the lattice structure, to confirm 
its validity further. 

It may be worth emphasizing once again that en- 
tropy/probability, which is the central concept of equi- 
librium thermodynamics/statistical mechanics, seems to 
have been left out of the investigations on nonequilibrium 
systems and pattern formation. Indeed, they have been 
almost always based on deterministic equations closely 
connected with conservation laws£ Calculations of en- 
tropy for open driven systems imply treating those finite 
systems as subjects of statistical mechanics with consid- 
ering the boundary conditions explicitly. Those calcula- 
tions are expected to shed new light on nonequilibrium 
phenomena in general which are still mysterious. 
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